CIE4604 - Simulation and Visualisation

Assignment 3 - Time Series Analysis of Vegetation Index

Student Name(s): Pratyush Kumar
Student Number(s): 5359252

In this assignment we do a time series analysis of NDVI using the HANTS (Harmonic Analysis of Time Series) algorithm and visualize the results in Python and QGIS.

For a description of the assignment consult the CIE4604-A3-NDVI_HANTS.pdf document.

Part I - Compute and Visualize NDVI in QGIS

If you prefer, you may add the QGIS results of part I here. See CIE4604-A3-NDVI_HANTS.pdf a detailed description of the actions.

In your report (or here) include at least the following results:

• Part I.a

  • Layout with the two NDVI maps, side by side, using the same map scale, with titles, legend and scale bar. You have to choose an appropriate representation (colormap) for the NDVI data (your grade depends on this)
  • A brief explanation on the method you used to extracted the data for the area of interest, with answers to the questions
  • Your motivation for the chosen representation (colormap) for the NDVI data

• Part I.b

  • Layout with the categorized training dataset.
  • Screenshot of the excel file with the attribute table

Furthermore, include in the zip file accompanying the report,

• The excel file with attributes (including zonal statistics)
• The re-projected training dataset (you will also need this for Part II)

Answer

For extraction of the raster images by the extent defined by the shapefile, i did the following:

1. reproject the bounding area to UTM Zone 30N
2. used GDAL to crop out the extent of the images as raster using the reprojected vector layer
3. used this cropped extent to extract the Raster images using GDAL's clip raster by mask layer
4. After this the raster calculator was used to produce the NDVI.

The colormap used for this purpose, is a diverging colormap with shades of brown for values below 0 and shades of green for values above 0. This was chosen also in a color blind freindly version of brown and green so that the maps are visible to everyone as intended. Moreover, the color scheme was chosen such that the greener pixel is indicative of higher NDVI value.

The GDAL command to make raster masks were done as follows

The image below is of the layout formed using QGIS:

The image below is of the attribute table exported as an xlsx file:

Part II - NDVI time series analysis with HANTS

HANTS is a software using the Fourier Analysis in order to detect outliers (clouds) and reconstruct image time series.
The analysis is performed over Sector BXII, an agricultural area located close to Sevilla (Spain). Input data files for this script are:

• NDVI_BXII_time_series_2017_2018_S2A_SP.tif (Geotiff with NDVI time series)
• NDVI_BXII_time_series_images_list.csv (Csv-file with the date, day number and day of year for the images)

Optional (non-required) input file(s) are:

• Training_2017_UTM30N_WGS84_SP.shp (Training dataset)

The actual HANTS algorithm is implemented in the Python function hants.py, see also help(hants.py).

This script is a template that you have to modify in places to obtain the desired results. Enter your answers as markdown.

# Importing all the libraries required for this assignment
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colorbar as colorbar
import matplotlib.colors as colors
import matplotlib.animation as animation
import matplotlib as mpl
import pandas as pd
import rasterio
import shapefile 

from scipy.io import loadmat
from datetime import date, datetime
from hants import hants

Load NDVI for Sector BXII, Sevilla, Spain

The input files are:

• NDVI_BXII_time_series_2017_2018_S2A_SP.tif - geotiff with NDVI
• NDVI_BXII_time_series_images_list.csv (Csv-file with the date, day number and day of year for the images)

The code for reading these dataset was already developed in Exercise 4. The following variables are created

• NDVI - 3-D array with dimensions [time, lat, lon] , i.e. NDVI[k,:,:] contains the k'th NDVI image at time ts[k]
• ts, tsdoy, tsdatestring - 1-D arrays with respectively number of days since 12/04/2017 (ts), day of year (tsdoy) and datestrings (tsdatestring).

You should definitely take a look at each variable to see what it is about.

# set the filenames and filepath
path = './input_files/Hants/'
NDVIfile = 'NDVI_BXII_time_series_2017_2018_S2A_SP.tif'
Datefile = 'NDVI_BXII_time_series_images_list.csv'

# open GeoTIFF dataset as an object 
src = rasterio.open(path + NDVIfile)

# print meta data (selection)
for key in src.meta:
  print(key,': ',src.meta[key])

print('\nbounds',src.bounds)
print('wkt',src.crs.wkt)

# compute imshow extent
extent = [ src.bounds.left ,  src.bounds.right ,  src.bounds.bottom ,  src.bounds.top ]
print('extent',extent)

# Read the image data
NDVI = src.read()
print('\nNDVI shape',NDVI.shape)

driver :  GTiff
dtype :  float64
nodata :  None
width :  866
height :  800
count :  33
crs :  EPSG:32630
transform :  | 20.00, 0.00, 214920.00|
| 0.00,-20.00, 4106200.00|
| 0.00, 0.00, 1.00|

bounds BoundingBox(left=214920.0, bottom=4090200.0, right=232240.0, top=4106200.0)
wkt PROJCS["WGS 84 / UTM zone 30N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-3],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32630"]]
extent [214920.0, 232240.0, 4090200.0, 4106200.0]

NDVI shape (33, 800, 866)

# Read dates using Pandas
datetable = pd.read_csv(path + Datefile) 

# Preview the first 5 lines of the loaded dates 
print("The date table .head() look like this: \n",datetable.head())

# Preview the data types
print(datetable.dtypes)

# Convert to 1-D numpy arrays 
ts = np.array(datetable.Day) 
tsdatestrings = np.array(datetable.ISOdate)
tsdoy = np.array(datetable.Doy)
print('ts',ts)
print('tsdatestrings',tsdatestrings)
print('tsdoy',tsdoy)
datetable

The date table .head() look like this: 
    Band      ISOdate  Day  Doy
0     1   2017-04-12    1  102
1     2   2017-04-22   11  112
2     3   2017-05-02   21  122
3     4   2017-05-12   31  132
4     5   2017-05-22   41  142
Band        int64
ISOdate    object
Day         int64
Doy         int64
dtype: object
ts [  1  11  21  31  41  51  61  71  81  91 101 111 131 141 151 161 171 181
 191 201 211 221 231 241 261 271 291 301 311 321 341 351 371]
tsdatestrings [' 2017-04-12' ' 2017-04-22' ' 2017-05-02' ' 2017-05-12' ' 2017-05-22'
 ' 2017-06-01' ' 2017-06-11' ' 2017-06-21' ' 2017-07-01' ' 2017-07-11'
 ' 2017-07-21' ' 2017-07-31' ' 2017-08-20' ' 2017-08-30' ' 2017-09-09'
 ' 2017-09-19' ' 2017-09-29' ' 2017-10-09' ' 2017-10-19' ' 2017-10-29'
 ' 2017-11-08' ' 2017-11-18' ' 2017-11-28' ' 2017-12-08' ' 2017-12-28'
 ' 2018-01-07' ' 2018-01-27' ' 2018-02-06' ' 2018-02-16' ' 2018-02-26'
 ' 2018-03-18' ' 2018-03-28' ' 2018-04-17']
tsdoy [102 112 122 132 142 152 162 172 182 192 202 212 232 242 252 262 272 282
 292 302 312 322 332 342 362   7  27  37  47  57  77  87 107]

	Band	ISOdate	Day	Doy
0	1	2017-04-12	1	102
1	2	2017-04-22	11	112
2	3	2017-05-02	21	122
3	4	2017-05-12	31	132
4	5	2017-05-22	41	142
5	6	2017-06-01	51	152
6	7	2017-06-11	61	162
7	8	2017-06-21	71	172
8	9	2017-07-01	81	182
9	10	2017-07-11	91	192
10	11	2017-07-21	101	202
11	12	2017-07-31	111	212
12	13	2017-08-20	131	232
13	14	2017-08-30	141	242
14	15	2017-09-09	151	252
15	16	2017-09-19	161	262
16	17	2017-09-29	171	272
17	18	2017-10-09	181	282
18	19	2017-10-19	191	292
19	20	2017-10-29	201	302
20	21	2017-11-08	211	312
21	22	2017-11-18	221	322
22	23	2017-11-28	231	332
23	24	2017-12-08	241	342
24	25	2017-12-28	261	362
25	26	2018-01-07	271	7
26	27	2018-01-27	291	27
27	28	2018-02-06	301	37
28	29	2018-02-16	311	47
29	30	2018-02-26	321	57
30	31	2018-03-18	341	77
31	32	2018-03-28	351	87
32	33	2018-04-17	371	107

Colormap selection

To plot the NDVI for the first day we are going to use matplotlib's imshow method to make a pseudo-color map. A good colormap for the NDVI is essential! It is up to you to come up with a good colormap for the NDVI, taking into account the lessons learned from Exercise 1 and 2. You may also use your answer from Exercise 4 here (this is a repetition).

Start by formulating the requirements for the color map. What kind of color map should we select? What is the range? Do we use a linear map, is it continuous or discrete, do we use two tones? What kind of normalization are we going to use? In any case, the default color map is certainly not a good choice. Why?

Have a look for instance at the following article https://publiclab.org/notes/cfastie/08-26-2014/new-ndvi-colormap to get some inspiration.

Explain here your considerations for the colormap and the reasons for choice in markdown

In the section below set the cmap variable with the chosen colormap.

In Python you can read the text file or mat file and create the colormap using matplotlib.colors.ListedColormap. See also Exercise 1 and 2.

cmapTable = scipy.io.loadmat('file.mat')['colors']
cmap = colors.ListedColormap(data/255)

You can then declare this new colormap in the imshow function by declaring the input as cmap = cmap

# SELECT APPROPRIATE COLORMAP FOR NDVI
cmap = 'RdYlGn'    # default colormap for matplotlib, change this for something better!!
# END COLORMAP SELECTION
clim = [-1 , 1]     # NDVI range

Pseudo color image plot of unfiltered NDVI for the first day

To plot the NDVI for the first day we are going to use matplotlib's imshow method to make a pseudo-color map. Use the colormap that you selected above using the cmap = cmap optional parameter in imshow, together with 'clim' and 'extent'.

The optional extent parameter is used to set the proper x- and y-axis coordinates. It is computed from the image bounds, while converting meter units to km.

# Put your code here ((you may borrow it from Exercise 4)

extentKiloMeter = np.array(extent)/1000
# borrowed from exercise 4
%matplotlib inline

band = 0

plt.figure(figsize=(8, 8))
fig = plt.imshow(NDVI[band, :, :], extent=extentKiloMeter, cmap=cmap, clim=clim)
plt.title('NDVI 2017-04-12')
plt.title('NDVI ' + tsdatestrings[band] + ' DAY ' + str(int(ts[band])))
plt.ylabel('North [km]')
plt.xlabel('East [km]')
plt.colorbar(label='NDVI [-]', fraction=0.0417, pad=0.06)
plt.show()

Create a video with the unfiltered NDVI data

A video is created showing the (unfiltered) NDVI data. We use imshow to display matrix NDVI(i,:,:) as an image for day i, using the colormap of your choice (see previous section).

The figure statement is called only once, imshow is called from within the loop over days, and the output is appended to the list of frames.

There are two ways to display the animation from a notebook,

• show the animation in a separate window using qt
• show the animation inline usign html5

You can also save the animation as either a mpeg file or an animated gif (both are commented out at the moment). The code is borrowed from Exercise 4.

# Put your code here, you may borrow it from Exercise 4
# // as is copy paste from exercise 4

# %matplotlib qt        # for using Qt
# %matplotlib notebook  # for Qt inline in notebook
# %matplotlib inline    # for inline html5

fig, ax=plt.subplots(figsize=(10,8))

# frames is a list of lists, each row is a list of artists to draw in the current frame; 
# here we are just animating one artist, the image, in each frame
frames = []
ndates = NDVI.shape[0]
for i in range(ndates):
    titlestr = 'NDVI ' + tsdatestrings[i] + ' DAY ' + str(int(ts[i]))   #     title = Add your title here
    im = ax.imshow(NDVI[i, :, :], extent=extentKiloMeter, aspect='equal', cmap=cmap, clim=clim, animated=True) #     im = Add your images here
    title = ax.text(0.5, 1.05, titlestr, size=plt.rcParams["axes.titlesize"], ha="center", transform=ax.transAxes, animated=True)
    frames.append([im, title])

plt.ylabel('North [km]')
plt.xlabel('East [km]')
fig.colorbar(im, ax=ax, label='NDVI [-]', fraction=0.0417, pad=0.06)

ani = animation.ArtistAnimation(fig, frames, interval=300, blit=False, repeat_delay=500)

FFwriter = animation.FFMpegWriter(fps=3) 
ani.save('NDVI_raw.mp4', writer = FFwriter)

The colormap used in the above animation and the plot is a diverging colormap from red to green. It has been chosen to portray that the greener the pixel is, the more NDVI value it has and vice versa with respect to the reds in the map.

Compute reconstructed time series using HANTS

Inputs for HANTS are the times ts (days) and the NDVI time series for one pixel. For the other input parameters see the help that comes with HANTS.

help(hants)

Help on function hants in module hants:

hants(ni, nb, nf, y, ts, HiLo, low, high, fet, dod, delta)
    Harmonic ANalysis of Time Series (HANTS)
    
    Parameters
    ----------
    ni :  int
        Number of images (total number of actual samples of the time series)
    nb : int 
        Length of the base period, measured in virtual samples (days, dekads, months, etc.)
    nf: int
        Number of frequencies to be considered above the zero frequency
    y: ndarray
        1D array of size size `ni` with input sample values (e.g. NDVI values)
    ts: ndarray
        1D array of size `ni` of time sample indicators (indicates virtual 
        sample number relative to the base period); numbers in array `ts` 
        maybe greater than `nb`
    HiLo: str {'Hi', 'Lo'}
        2-character string indicating rejection of high 'Hi' or low 'Lo' outliers
    low: float
        valid range minimum
    high: float
        valid range maximum (values outside the valid range are rejeced right away)
    fet: float
        fit error tolerance (points deviating more than fet from curve fit are rejected)
    dod : int
        degree of overdeterminedness (iteration stops if number of points reaches 
        the minimum required for curve fitting, plus dod). This is a safety measure
    delta : float
        small positive number (e.g. 0.1) to suppress high amplitudes
    
    Returns
    -------
    amp : ndarray
        1D array of length `nf`+` with Amplitude for frequencies 0...nf
    phi : ndarray
        1D array of length `nf'+1 with Phase in degrees for frequencies 0...nf
    yr : ndarray
        1D array of length `ni` with smoothed values (e.g. NDVI)   
    outliers : ndarray 
        1D Boolean array of length `ni` with flagged outliers
     
    
    Notes
    -----    
    This function applies the Harmonic ANalysis of Time Series (HANTS)
    algorithm originally developed by the Netherlands Aerospace Centre (NLR)
    (http://www.nlr.org/space/earth-observation/).
    
    This python implementation was based on two previous implementations
    available at the following links:
    https://codereview.stackexchange.com/questions/71489/harmonic-analysis-of-time-series-applied-to-arrays
    http://nl.mathworks.com/matlabcentral/fileexchange/38841-matlab-implementation-of-harmonic-analysis-of-time-series--hants-
    
    Original Author for Python Version: Espinoza-Dávalos, G. E., Bastiaanssen, W. G. M., Bett, B., & Cai, X. (2017).
    A Python Implementation of the Harmonic ANalysis of Time Series (HANTS) Algorithm for Geospatial Data.
    http://doi.org/10.5281/zenodo.820623
    
    Edited by: Ullas Rajvanshi and Hans van der Marel

Based on the help information define your HANTS input parameters below. Some are preset, others are left open.

nb= 371       # length of the base period (usually one year) measured in timeunits of ts (number of days)
nf= 1         # number of frequencies to be considered above the zero frequency ( 1 ... 6)
fet= 0.05      # fit error tolerance (points deviating more than fet from curve fit are rejected, 0.05 ... 0.3)
HiLo='Lo'     # 2-character string indicating rejection of high or low outliers ('Lo' or 'Hi')
low=-1        # valid range minimum (values below are rejected right away)
high=+1       # valid range maximum (values above are rejected right away)
dod=1         # degree of overdeterminedness (iteration stops if number of points reaches the minimum required)
delta=0.1     # small positive number (e.g. 0.1) to suppress high amplitudes

Set the line and sample number of the pixel for which to apply HANTS.

We will show you later how to find the line and sample numbers for specific pixel areas, e.g. one of the field in the training set.

row = 348
col = 299

Do HANTS, first we have to extract the pixel data into the array datain. We use the Matlab function squeeze to do this, because NDVI is stored as a three dimensional array and we have to collect all values that belong to a single pixel. We use squeeze because NDVI[:,row,col] does not return an array, but a list, and squeeze will convert this list to an array.

HANTS Outputs are:

  amp   = array of amplitudes, first element is the average of the curve
  phi   = array of phases, first element is zero
  dataout   = array holding reconstructed time series 
  outliers  = boolean array with the outliers 

ndates = NDVI.shape[0]
datain = np.squeeze(NDVI[:, row, col])

amp, phi, dataout, outliers = hants(ndates, nb, nf, datain, ts, HiLo, low, high, fet, dod, delta)

print('amp',amp)
print('phi',phi)
print('dataout',dataout)
print('outliers',outliers)
numoutliers=np.sum(outliers)
print('numoutliers',numoutliers)

amp [0.49355339 0.40122089]
phi [ 0.         51.22445148]
dataout [0.75008759 0.79841392 0.83801707 0.86776384 0.88680308 0.89459001
 0.8909018  0.87584399 0.84984744 0.81365601 0.76830527 0.71509287
 0.59135488 0.52436988 0.45650312 0.3896965  0.3258616  0.26682499
 0.2142759  0.16971797 0.13442617 0.10941031 0.09538619 0.09275511
 0.121645   0.15233934 0.24186053 0.29812585 0.35998307 0.42566224
 0.5609135  0.62661556 0.74482646]
outliers [ True  True  True  True  True  True False False False False  True  True
  True  True  True  True  True  True  True  True False False  True False
  True  True False  True  True  True  True  True  True]
numoutliers 25

Plot the original vs reconstructed time series

%matplotlib inline


# making a function to replicate the plotting
def readymade_HANTS_plot_maker(ts, datain, outliers, dataout, row, col, nb, nf, fet):
    numoutliers=np.sum(outliers)
    plt.figure(figsize=(10,6))
    plt.plot(ts, datain, 'b*', label='Original NDVI')
    plt.plot(ts[outliers], datain[outliers], 'rx', markersize=10, markeredgewidth=3, label='Outlier NDVI')
    plt.plot(ts, dataout, 'r', label='Smoothed with HANTS')
    plt.xlim([1, 365])
    #plt.ylim([0, 1])
    plt.xlabel('Days from 12 April 2017')
    plt.ylabel('NDVI [-]')
    plt.title('NDVI - Sentinel 2 (row {}, column {})'.format(row, col));
    plt.text(6,0.77,'nb={}\nnf={}\nfet={}\noutliers={}'.format(nb, nf, fet, numoutliers))
    plt.legend()
    plt.show()

readymade_HANTS_plot_maker(ts, datain, outliers, dataout, row, col, nb, nf, fet)

Plot HANTS for different HANTS params

this will help us to iterate through different values of frequencies and days

# check hants with different frequencies
def HANTS_and_plottting(ndates, nb, nf, datain, ts, HiLo, low, high, fet, dod, delta, row,col):
    ndates = NDVI.shape[0]
    datain = np.squeeze(NDVI[:, row, col])

    amp, phi, dataout, outliers = hants(ndates, nb, nf, datain, ts, HiLo, low, high, fet, dod, delta)
    numoutliers=np.sum(outliers)
    readymade_HANTS_plot_maker(ts, datain, outliers, dataout, row, col, nb, nf, fet)
    

Plotting for different vegetation pixels

HANTS_and_plottting(ndates, nb, 3, datain, ts, HiLo, low, high, 10/100, dod, delta, row, col)
HANTS_and_plottting(ndates, nb, 3, datain, ts, HiLo, low, high, 10/100, dod, delta, 200,150)

# for i in range(1, 7, 1):
#     # changing frequency for everything remaining same
#     for j in range(5, 31, 5) :
#         HANTS_and_plottting(ndates, nb, i, datain, ts, HiLo, low, high, j/100, dod, delta)
#     print(f"Changing frequency to {i+1}\n")
#     # HANTS_and_plottting(ndates, nb, nf, datain, ts, HiLo, low, high, , dod, delta)a

Select suitable HANTS input parameters

Use the above code as a template to select a suitable set of input parameters. You can try different values for the number of frequencies (nf) , base period (nb) and fit error tolerance.

You should at least

• Plot the reconstructed vs the original time series for two vegetated pixels characterized by different temporal trend (different crop growing cycles).
• Plot the reconstructed vs the original time series for different HANTS input parameters, but ONLY GIVE RESULTS FOR A FEW SENSIBLE SELECTIONS to underline your choice for the frequency and/or error tolerance.

In the report, explain your choice for the base period (you don't have to show plots, there is only one good value for the base period), and include the plots for the chosen parameters.

# Add your plots here
nb=371
nf=3
fet=0.1

Base period: It is chosen as 371 since the 33 NDVI images are over the period of 371 days.
nf: 3
fet: 0.1
From the above plots, it is decided to take nf value as 3 and the tolerance set at 0.1, since at nf=4 we start seeing that the curve almost fits the highest values whereas on increasing the nf value further it is observed that the curve can be called as over-fitting the values.

Apply HANTS for all pixels

Apply the HANTS algorithm for all pixels, and save the amplitude, phase and reconstucted NDVI.

THIS IS A VERY TIME CONSUMING PROCESS THAT MAY TAKE 15-30 MINUTES! SO ONLY DO THIS WHEN YOUR ARE SURE YOU GOT THE CORRECT PARAMETERS, OR TEST IT FIRST ON A SMALLER SAMPLE.

print('Smoothing and filling the missing data, please be patient ... ')


# Initialize new arrays to collect the HANTS output
ni, ny, nx = NDVI.shape
NDVI_HANTS = np.zeros((ni, ny, nx))
amp = np.zeros((nf + 1, ny, nx))
phi = np.zeros((nf + 1, ny, nx))

# Call HANTS in a double loop, line by line (SEE THE WARNING)
#for row in range(0, 10): # use this for smaller sample
for row in range(0, ny):
    if ( row % 25 ) == 0: print("Doing line {} out of {}".format(row, ny))  
    for col in range(0, nx):
        data = NDVI[:, row, col]
        if np.isnan(data).sum() != ni:
            data[np.isnan(data)] = low - 1.0
            amp_, phi_, yr_, outlier_ = hants(ndates, nb, nf, data, ts, HiLo, low, high, fet, dod, delta)
            amp[:, row, col], phi[:, row, col], NDVI_HANTS[:, row, col] = amp_, phi_, yr_

print("Done with smoothing and filling")

Smoothing and filling the missing data, please be patient ... 
Doing line 0 out of 800
Doing line 25 out of 800
Doing line 50 out of 800
Doing line 75 out of 800
Doing line 100 out of 800
Doing line 125 out of 800
Doing line 150 out of 800
Doing line 175 out of 800
Doing line 200 out of 800
Doing line 225 out of 800
Doing line 250 out of 800
Doing line 275 out of 800
Doing line 300 out of 800
Doing line 325 out of 800
Doing line 350 out of 800
Doing line 375 out of 800
Doing line 400 out of 800
Doing line 425 out of 800
Doing line 450 out of 800
Doing line 475 out of 800
Doing line 500 out of 800
Doing line 525 out of 800
Doing line 550 out of 800
Doing line 575 out of 800
Doing line 600 out of 800
Doing line 625 out of 800
Doing line 650 out of 800
Doing line 675 out of 800
Doing line 700 out of 800
Doing line 725 out of 800
Doing line 750 out of 800
Doing line 775 out of 800
Done with smoothing and filling

Save the results now!!

Variables can be saved and loaded with the Spyder IDE. Do that for the newly computed results, so that for the follow on work you don't have to redo the time consuming processing every time. Here we will make use of pickle library to store the variables so that we can always reload them later.

# # Save important variables
import pickle

# f = open('NDVIr.pckl', 'wb')
# pickle.dump([ts, tsdatestrings, extent, NDVI_HANTS, amp, phi], f)
# f.close()

If the variables have to be reloaded you can use the following code.

# ReLoad the variables
f = open('NDVIr.pckl', 'rb')
ts, tsdatestrings, extent, NDVI_HANTS, amp, phi = pickle.load(f)
f.close()

Write GeoTIFF files with reconstructed NDVI, amplitude and phase

Three new GeoTIFF files are created, with

• the reconstructed (smoothed) NDVI
• the amplitude
• the phase

These files can easily be imported by other software such as QGIS.

# Save reconstructed NDVI as GeoTiff file
profile = src.profile
with rasterio.open('NDVI_reconstructed.tif', 'w', **profile) as dst:
    for id, layer in enumerate(NDVI_HANTS, 1):
        dst.write_band(id, NDVI_HANTS[id - 1, :, :])

# Save Amplitude as GeoTiff file
profile = src.profile
profile['count']=len(amp)
with rasterio.open('NDVI_amplitude.tif', 'w', **profile) as dst:
    for id, layer in enumerate(amp, 1):
        dst.write_band(id, amp[id - 1, :, :])

# Save Phase as GeoTiff file
profile = src.profile
profile['count']=len(phi)
with rasterio.open('NDVI_phase.tif', 'w', **profile) as dst:
    for id, layer in enumerate(phi, 1):
        dst.write_band(id, phi[id - 1, :, :])
        

Create a video of the reconstruced NDVI data

Create a video showing the reconstructed NDVI data.

# Your code here

# create a function to return a list of image and title
def dayMapper(no_day):
    reconst_date_ = tsdatestrings[no_day]
    img_local_ax = ax.imshow( NDVI_HANTS[no_day, :, :], cmap=cmap, clim=clim, animated=True, aspect="equal" )
    img_title = ax.text(150, -20, f"NDVI Reconstructed with HANTS on {reconst_date_}")
    return [img_local_ax, img_title]


fig, ax = plt.subplots()
ndates = NDVI_HANTS.shape[0]

# make frames
frames = [dayMapper(i) for i in range(1, ndates)]

# generate animation object
anim = animation.ArtistAnimation(    fig, frames, interval=300, blit=False, repeat_delay=500 )

# save resultsF
fig.colorbar(frames[0][0], ax=ax, fraction=0.0417, pad=0.06)
fig.set_facecolor("white")


anim.save('NDVI_reconstructed_animation.mp4', writer = FFwriter)

Instead of an animation we can also create a lot of subplots using the code below.

%matplotlib inline
ndates = NDVI_HANTS.shape[0]

ncols = 5
nrows = int(np.ceil(ndates/ncols))

fig, axs = plt.subplots(nrows, ncols, figsize=(21,27), tight_layout=True)
axs= axs.flat
[axi.set_axis_off() for axi in axs]

for i in range(ndates):
    titlestr = 'Reconstructed NDVI ' + tsdatestrings[i] + ' DAY ' + str(int(ts[i]))
    axs[i].imshow(NDVI_HANTS[i, :, :], extent=extent, aspect='equal', cmap=cmap, clim=clim)
    axs[i].set_xticks([])
    axs[i].set_yticks([])
    axs[i].set_title(titlestr)

plt.show()

Plot amplitude and phase from HANTS (Part 2.D)

Make a figure with in the subplots a pseudo color image of the amplitude and phase for each frequency. Please think carefully about the colormap to use.

Which colormap do you use for the amplitude. Should the colormap for the phase be the same, or different. Please explain why.

# Your code here

plt.figure(figsize=(15,15))
# first plot
plt.subplot(1, 2, 1)
amplitude_cmap = 'PuBuGn'#'YlGnBu'
plt.imshow(amp[1, :, :], cmap=amplitude_cmap, label="Amplitude", extent=extentKiloMeter)
plt.colorbar(fraction=0.042)
plt.ylabel("North [km]")
plt.xlabel("East [km]")
plt.title("Amplitude of reconstructed time series")

# 2nd plot
plt.subplot(1, 2, 2)
phase_cmap = 'YlGnBu'#'PuBuGn'
plt.imshow(phi[1, :, :], cmap=phase_cmap, extent=extentKiloMeter)
plt.colorbar(fraction=0.042, label="Phase angle [Degrees]")
plt.title("Phase angle of reconstructed time series")
plt.ylabel("North [km]")
plt.xlabel("East [km]")

plt.show()

Plotting cloudy vs reconstructed NDVI from a cloudy day (Part 2.E)

28 November 2017 was significantly cloudy, so we will use that image, from the dataframe we can determine that the index for this image is 22

plt.figure(figsize=(15,15))
# first plot
plt.subplot(1, 2, 1)
plt.imshow(NDVI[22, :, :], cmap=cmap, extent=extentKiloMeter)
plt.colorbar(fraction=0.042)
plt.ylabel("North [km]")
plt.xlabel("East [km]")
plt.title("Unfiltered NDVI 2017-11-28")

# 2nd plot
plt.subplot(1, 2, 2)
plt.imshow(NDVI_HANTS[22, :, :], cmap=cmap, extent=extentKiloMeter)
plt.colorbar(fraction=0.042, label="NDVI")
plt.ylabel("North [km]")
plt.xlabel("East [km]")
plt.title("Reconstructed NDVI 2017-11-28")

plt.show()

Analyze the training dataset

The training dataset is a shapefile with information on the actual crop that is growing in the fields. The original training dataset was in UTM zone 29N. We used QGIS to convert this to UTM zone 30N, which was is also the zone that is used for the NDVI data.

First we must read the shape file with training data

training = shapefile.Reader("Training_set/Training_2017_UTM30N_WGS84_SP.shp")

# Your code here

training = shapefile.Reader("./input_files/Training_set/Training_2017_UTM30N_WGS84_SP.shp")
training

<shapefile.Reader at 0x7f46f4cefa60>

The training dataset consists of 15 polygons. We will create an plot with the amplitude data and overlay the polygon boundaries.

plt.figure(figsize=(12,8))
for shape in training.shapeRecords():
    x = [i[0]/1000 for i in shape.shape.points[:]]
    y = [i[1]/1000 for i in shape.shape.points[:]]
    plt.plot(x,y, color='yellow', linewidth=2)
greys = mpl.cm.get_cmap('Greys_r')
plt.imshow(amp[0, :, :], extent=extentKiloMeter, cmap=greys)
plt.axis('image')
plt.colorbar(label='Amplitude [-]')
plt.ylabel('North [km]')
plt.xlabel('East [km]')
plt.title('NDVI Training set')
plt.show()

Part III - Visualize HANTS output in QGIS

If you prefer, you may add the QGIS results of part III here. See CIE4604-A3-NDVI_HANTS.pdf a detailed description of the actions.

The amplitude obtained from the reconstruction of the images was loaded into QGIS and the amplitudes stored in bands 1,3,4 were mapped to the R,G,B bands on the visible image.

The image finally looks like the following:

[End of the assignment]